﻿// 标准库或者第三方库

// 自定义头文件
#include "template_problem01.hxx"

// ACIS
#include "acis/acis.hxx"
#include "acis/api.err"
#include "acis/api.hxx"
#include "acis/body.hxx"
#include "acis/check.hxx"
#include "acis/cstrapi.hxx"
#include "acis/logical.h"
#include "acis/module.hxx"
#include "acis/primtive.hxx"

#include "acis/law_base.hxx"
#include "acis/law.hxx"
#include "acis_utils.hpp"
#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>

using namespace boost::multiprecision;


// 自定义比较函数，考虑精度问题
bool are_equal_with_precision(cpp_dec_float_50 a, cpp_dec_float_50 b, cpp_dec_float_50 epsilon = 1e-9) {
	return fabs(a - b) < epsilon;
}

// 自定义 unique 比较函数
struct PrecisionComparator {
	cpp_dec_float_50 epsilon;
	PrecisionComparator(cpp_dec_float_50 eps = 1e-9) : epsilon(eps) {}

	bool operator()(cpp_dec_float_50 a, cpp_dec_float_50 b) const {
		return are_equal_with_precision(a, b, epsilon);
	}
};

/**
 * @brief	判断一个数是否在给定精度下接近0
 * @note	使用高精度库
 * @param[in]	number	待检测的数字
 * @param[in]	tol		给定的精度
 * @param[out]			是否接近0
 */
logical near_zero_2(cpp_dec_float_50 number, cpp_dec_float_50 tol = 1e-30)
{
	return fabs(number) <= tol;
}
/**
 * @brief		二次函数求根
 * @param[in]	f		二次函数表达式
 * @return		roots	二次函数的根(0个 / 1个 / 2个)
 */
std::vector<double> twice_order_roots(law* f);

/**
 * @brief		二次函数求根
 * @param[in]	coef	二次函数系数数组
 * @param[in]	N		系数个数
 * @return		roots	二次函数的根(0个 / 1个 / 2个)
 */
std::vector<cpp_dec_float_50> twice_order_roots_2(cpp_dec_float_50 coef[], int N);

/**
 * @brief		牛顿切线法迭代公式
 * @detail		x_next = x - f(x)/f'(x)
 * @param[in]	f_x		目标函数值
 * @param[in]	df_x	目标函数的导数值
 * @param[in]	x		迭代初始点
 * @return		x_next	迭代下一个点
 */
inline double newton_iter(law* f, law* df, double x);

/**
 * @brief		三阶改进牛顿法迭代公式
 * @param[in]	f_x		目标函数值
 * @param[in]	df_x	目标函数的一阶导数值
 * @param[in]	ddf_x	目标函数的二阶导数值
 * @param[in]	x		迭代初始点
 * @param[in]	m		权重因子
 * @return		x_next	迭代下一个点
 */
inline double three_order_newton(double f_x, double df_x, double ddf_x, double x, double m = 15.0);

/**
 * @brief 判断函数在一个单调区间内是否有零点
 * @param[in]	f			函数表达式
 * @param[in]	start		区间左端点
 * @param[in]	end			区间右端点
 * @return		true/false	是否存在零点
 */
bool hasroot(law* f, double start, double end);

/**
 * @brief	判断函数在一个单调区间内是否有零点
 * @note	使用高精度库
 * @param[in]	coef		多项式函数系数数组
 * @param[in]	start		区间左端点
 * @param[in]	end			区间右端点
 * @param[in]	N			多项式系数个数
 * @return		true/false	是否存在零点
 */
bool hasroot_2(cpp_dec_float_50 coef[], cpp_dec_float_50 start, cpp_dec_float_50 end, int N);

/**
 * @brief	在单调区间内求根(求根算法的核心)
 * @detail	使用牛顿切线法求根，超出边界时使用牛顿割线法拉回区间
 * @param[in]	f		函数表达式
 * @param[in]	start	区间左端点
 * @param[in]	end		区间右端点
 * @return		root	零点
 *
 * @note [start，end]区间应该是单调的，需要先使用hasroot确保在区间内有根，仅存在一个零点
 */
double getroot(law* f, double start, double end);

/**
 * @brief	求出多项式方程在区间内的所有零点
 * @detail	将多项式函数划分为单调区间，在每一个单调区间内部求解单根
 * @param[in]	f		多项式函数
 * @param[in]	start	区间左端点
 * @param[in]	end		区间右端点
 * @return		返回区间内所有的根
 */
std::vector<double> getroots(law* f, double start, double end);

/**
 * @brief	求出多项式方程的单调区间
 * @detail	求导函数的零点，并且划分导函数的正区间、负区间，即为原函数的增减区间
 * @return	存放导函数零点的数组
 *			[start, x1, x2, x3, x4, ... , xn, end]
 *			任意相邻点构成的区间内，原函数为增函数/减函数
 *			e.g. [start, x1]减区间，[x_i, x_(i+1)]增区间 ...
 */
std::vector<double> get_mon_interval(law* f, double start, double end);

/**
 * @brief	系数多项式代入x求值
 * @note	使用高精度库
 * @param[in]	coef	多项式系数数组
 * @param[in]	x		代入的函数点x
 * @param[in]	N		系数个数
 * @param[out]			函数值
 */
cpp_dec_float_50 polynomial_eval(cpp_dec_float_50 coef[], cpp_dec_float_50 x, int N);

/**
 * @brief	系数多项式求导
 * @note	使用高精度库
 * @param[out]	deriv	导数多项式系数数组
 * @param[in]	coef	原多项式系数数组
 * @param[in]	N		系数个数
 */
void polynomial_derivative(cpp_dec_float_50 deriv[], const cpp_dec_float_50 coef[], int N);

/**
 * @brief	求出多项式方程的单调区间
 * @note	使用高精度库
 * @detail	求导函数的零点，并且划分导函数的正区间、负区间，即为原函数的增减区间
 * @return	存放导函数零点的数组
 *			[start, x1, x2, x3, x4, ... , xn, end]
 *			任意相邻点构成的区间内，原函数为增函数/减函数
 *			e.g. [start, x1]减区间，[x_i, x_(i+1)]增区间 ...
 */
std::vector<cpp_dec_float_50> get_mon_interval_2(cpp_dec_float_50 coef[], cpp_dec_float_50 start, cpp_dec_float_50 end, int N);

/**
 * @brief	求出多项式方程在区间内的所有零点
 * @detail	将多项式函数划分为单调区间，在每一个单调区间内部求解单根
 * @note	使用高精度库
 * @param[in]	coef	多项式函数系数数组
 * @param[in]	start	区间左端点
 * @param[in]	end		区间右端点
 * @param[in]	N		系数个数
 * @return		返回区间内所有的根
 */
std::vector<cpp_dec_float_50> getroots_2(cpp_dec_float_50 coef[], cpp_dec_float_50 start, cpp_dec_float_50 end, int N);


std::vector<double> twice_order_roots(law* f) {
	std::vector<double> roots = {};

	API_BEGIN
		law_polynomial* f_poly = f->polynomial(f);
	double a = f_poly->coefficent(2);
	double b = f_poly->coefficent(1);
	double c = f_poly->coefficent(0);
	// 确保是二次多项式
	if (f_poly->degree() != 2) sys_error(API_FAILED);

	// 使用求根公式
	double delta = pow(b, 2) - (4 * a * c);
	if (delta > 0) {
		roots.push_back((-b - sqrt(delta)) / (2 * a));
		roots.push_back((-b + sqrt(delta)) / (2 * a));
	}
	else if (delta == 0) {
		roots.push_back(-b / (2 * a));
	}
	API_END

		return roots;
}



std::vector<cpp_dec_float_50> twice_order_roots_2(cpp_dec_float_50 coef[], int N) {
	std::vector<cpp_dec_float_50> roots = {};

	API_BEGIN

		cpp_dec_float_50 a = coef[2];
	cpp_dec_float_50 b = coef[1];
	cpp_dec_float_50 c = coef[0];
	// 确保是二次多项式
	if (N != 2) sys_error(API_FAILED);

	// 使用求根公式
	cpp_dec_float_50 delta = pow(b, 2) - (4 * a * c);
	//std::cout << "delta:" << delta << std::endl;
	if (delta > 0) {
		roots.push_back((-b - sqrt(delta)) / (2 * a));
		roots.push_back((-b + sqrt(delta)) / (2 * a));
	}
	else if (near_zero_2(delta)) {
		roots.push_back(-b / (2 * a));
	}
	API_END

		return roots;
}


cpp_dec_float_50 polynomial_eval(cpp_dec_float_50 coef[], cpp_dec_float_50 x, int N) {
	// 最高次项的系数
	cpp_dec_float_50 r = coef[N];
	// 使用Horner法则计算在X处的值
	for (int i = N - 1; i >= 0; --i) r = r * x + coef[i];
	return r;
}


void polynomial_derivative(cpp_dec_float_50 deriv[], const cpp_dec_float_50 coef[], int N) {
	deriv[0] = coef[1];
	for (int i = 2; i <= N; ++i) deriv[i - 1] = i * coef[i];
}

inline double newton_iter(double f_x, double df_x, double x) {
	if (near_zero_2(df_x)) {
		std::cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<< Warnning:Devide by zero <<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl;
		return x;
	}
	return x - f_x / df_x;;
}

inline cpp_dec_float_50 newton_iter_2(cpp_dec_float_50 f_x, cpp_dec_float_50 df_x, cpp_dec_float_50 x) {
	if (near_zero_2(df_x)) {
		std::cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<< Warnning:Devide by zero <<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl;
		return x;
	}
	return x - f_x / df_x;;
}

//double high_order_newton(law* f, double x) {
//
//	law* df = f->derivative(0);
//	law* d2f = df->derivative(0);
//	law* d3f = d2f->derivative(0);
//	law* d4f = d3f->derivative(0);
//	law* d5f = d4f->derivative(0);
//	law* d6f = d5f->derivative(0);
//	law* d7f = d6f->derivative(0);
//
//	double f_x = f->eval(x);
//	double df_x = df->eval(x);
//	double d2f_x = d2f->eval(x);
//	double d3f_x = d3f->eval(x);
//	double d4f_x = d4f->eval(x);
//	double d5f_x = d5f->eval(x);
//	double d6f_x = d6f->eval(x);
//	double d7f_x = d7f->eval(x);
//
//	double x_next;
//	double epsilon = 1e-6;
//	double eta = 0.6;
//
//
//	if (abs(f_x) >= epsilon) {
//		if (df_x != 0) {
//			x_next = x - f_x / df_x;
//		}
//		else {
//			x_next = x - eta * f_x;
//		}
//	}
//	else if (abs(df_x) >= epsilon) {
//		if (d2f_x != 0) {
//			x_next = x - 2 * df_x / d2f_x;
//		}
//		else {
//			x_next = x - eta * df_x;
//		}
//	}
//	else if (abs(d2f_x) >= epsilon) {
//		if (d3f_x != 0) {
//			x_next = x - 3 * d2f_x / d3f_x;
//		}
//		else {
//			x_next = x - eta * d2f_x;
//		}
//	}
//	else if (abs(d3f_x) >= epsilon) {
//		if (d4f_x != 0) {
//			x_next = x - 4 * d3f_x / d4f_x;
//		}
//		else {
//			x_next = x - eta * d3f_x;
//		}
//	}
//	else if (abs(d4f_x) >= epsilon) {
//		if (d5f_x != 0) {
//			x_next = x - 5 * d4f_x / d5f_x;
//		}
//		else {
//			x_next = x - eta * d4f_x;
//		}
//	}
//	else if (abs(d5f_x) >= epsilon) {
//		if (d6f_x != 0) {
//			x_next = x - 6 * d5f_x / d6f_x;
//		}
//		else {
//			x_next = x - eta * d5f_x;
//		}
//	}
//	else if (abs(d6f_x) >= epsilon) {
//		if (d7f_x != 0) {
//			x_next = x - 7 * d6f_x / d7f_x;
//		}
//		else {
//			x_next = x - eta * d6f_x;
//		}
//	}
//
//	return x_next;
//}

inline double three_order_newton(double f_x, double df_x, double ddf_x, double x, double m) {
	double num = 2 * pow(m, 2) * pow(f_x, 2) * ddf_x;
	double den = m * (3 - m) * f_x * df_x * ddf_x + pow(m - 1, 2) * pow(df_x, 3);
	if (near_zero(den)) {
		std::cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<< Warnning:Devide by zero <<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl;
		return x;
	}
	return x - num / den;
}

//super-Halley法
inline double four_order_newton_halley(double f_x, double df_x, double ddf_x, double ddf_x_2, double x) {
	double num = 6 * pow(df_x, 3) * f_x - ddf_x_2 * pow(f_x, 2);
	double den = 6 * pow(df_x, 4) - 2 * ddf_x_2 * df_x * f_x;
	if (near_zero(den)) {
		std::cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<< Warnning:Devide by zero <<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl;
		return x;
	}
	return x - num / den;
}

// Jarratt 法
inline double four_order_newton_Jarratt(double f_x, double df_x, double df_y, double x) {
	double num = 3 * df_y * f_x + df_x * f_x;
	double den = 6 * df_y * df_x - 2 * df_x * df_x;
	if (near_zero(den)) {
		std::cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<< Warnning:Devide by zero <<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl;
		return x;
	}
	return x - num / den;
}

bool hasroot(law* f, double start, double end) {
	if (start > end) {
		std::cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<< Internal ERROR <<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl;
		return false;
	}
	bool flag = false;
	API_BEGIN
		double f_a = f->eval(start), f_b = f->eval(end);
	if (near_zero_2(f_a) || near_zero_2(f_b) || (f_a * f_b) < 0) {
		flag = true;
	}
	API_END

		return flag;
}


double getroot(law* f, double start, double end) {
	double root;

	//API_BEGIN
	if (start > end) {
		std::cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<< Internal ERROR <<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl;
		//sys_error(API_FAILED);
		return root;
	}

	// 根在区间端点，直接返回
	if (near_zero_2(f->eval(start))) {
		return root = start;
	}
	if (near_zero_2(f->eval(end))) {
		return root = end;
	}
	//API_END

	// 根在区间内部
	API_BEGIN
		// 获取多项式的阶数
		law_polynomial* f_poly = f->polynomial(f);
	int order = f_poly->degree();

	// 一阶多项式ax+b=0，根root = -b/a
	if (order == 1) {
		root = -f_poly->coefficent(0) / f_poly->coefficent(1);
	}

	// 二阶多项式ax^2+bx+c=0
	else if (order == 2) {
		//std::vector twice_roots = twice_order_roots(f);
		for (double x : twice_order_roots(f)) {
			if (start <= x && x <= end) {
				root = x;
				break;
			}
		}
	}

	// 高阶多项式
	else {
		// 迭代初始点设置为中点
		double x = (start + end) / 2;
		double x_next;

		// 使用牛顿迭代法求根
		law* df = f->derivative(0);
		law* ddf = df->derivative(0);
		for (int i = 0; i < 2000; i++) {
			double f_x = f->eval(x);
			double df_x = df->eval(x);
			double ddf_x = ddf->eval(x);
			double ddf_x_2 = ddf->eval(x - f_x);	//super halley

			// 函数值足够小视为零点
			if (near_zero_2(f_x)) {
				root = x;
				break;
			}
			//if (fabs(f_x) > 0.01) {
			//	x_next = newton_iter(f_x, df_x, x);
			//}
			//else {
			//	// 三阶改进牛顿法
			//	//x_next = three_order_newton(f_x, df_x, ddf_x, x, 15.0);

			//	//四阶牛顿法super halley
			//	x_next = four_order_newton_halley(f_x, df_x, ddf_x, ddf_x_2, x);

			//	//四阶牛顿法Jarratt
			//	//x_next = four_order_newton_Jarratt(f_x, df_x, df_y, x);
			//}
			x_next = newton_iter(f_x, df_x, x);
			//x_next = high_order_newton(f, x);
			//x_next = four_order_newton_halley(f_x, df_x, ddf_x, ddf_x_2, x);
			// 步长过小，迭代收敛
			if (near_zero_2(x_next - x)) {
				root = x_next;
				break;
			}
			//缩小区间
			double f_xmin = f->eval(start);
			double f_xmax = f->eval(end);
			double f_xnext = f->eval(x_next);
			if (f_xnext * f_xmin > 0) {
				start = x_next;
			}
			if (f_xnext * f_xmax > 0) {
				end = x_next;
			}
			// 超出求根区间，使用牛顿切线法拉回
			if (x_next < start || x_next > end) {
				double f_a = f->eval(start);
				double f_x = f->eval(x);
				double f_b = f->eval(end);
				// 根在迭代点左侧
				if ((f_a * f_x) < 0) {
					x_next = start - f_a * (x - start) / (f_x - f_a);
				}
				// 根在迭代点右侧
				else {
					x_next = end - f_b * (x - end) / (f_x - f_b);
				}
			}
			x = x_next;

		}
		root = x;		// 否则陷入循环的时候会返回未初始化的root值
	}
	API_END
		return root;
}




std::vector<double> get_mon_interval(law* f, double start, double end) {
	std::vector<double> intervals;
	intervals.push_back(start);

	// 递归将导数的零点放入
	std::vector<double> df_roots = getroots(f->derivative(0), start, end);
	intervals.insert(intervals.end(), df_roots.begin(), df_roots.end());

	intervals.push_back(end);
	return intervals;
}


std::vector<double> getroots(law* f, double start, double end) {
	std::vector<double> roots;

	API_BEGIN
		// 获取多项式的阶数
		law_polynomial* f_poly = f->polynomial(f);
	int order = f_poly->degree();

	// 简单的一阶和二阶多项式可以直接求出
	// 一阶多项式ax+b=0，根root = -b/a
	if (order == 1) {
		double x = -f_poly->coefficent(0) / f_poly->coefficent(1);
		if (start <= x && x <= end) {
			roots.push_back(x);
		}
	}

	// 二阶多项式ax^2+bx+c=0
	else if (order == 2) {
		//std::vector twice_roots = twice_order_roots(f);
		for (double x : twice_order_roots(f)) {
			if (start <= x && x <= end) {
				roots.push_back(x);
			}
		}
	}

	// 高阶多项式
	else {
		// 求单调区间
		std::vector<double> intervals = get_mon_interval(f, start, end);
		/*std::cout << "---------------------------------已分区间---------------------------------" << std::endl;

		for (int i = 0; i < intervals.size() - 1; i++) {
			std::cout << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10) << "(" << intervals[i] << "," << intervals[i + 1] << ")" << std::endl;
		}*/
		// 对每个区间求根
		for (int i = 0; i < intervals.size() - 1; i++) {
			if (hasroot(f, intervals[i], intervals[i + 1])) {
				roots.push_back(getroot(f, intervals[i], intervals[i + 1]));
			}
		}
	}

	API_END
		return roots;
}


bool hasroot_2(cpp_dec_float_50 coef[], cpp_dec_float_50 start, cpp_dec_float_50 end, int N) {
	if (start > end) {
		std::cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<< Internal ERROR <<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl;
		return false;
	}
	bool flag = false;
	API_BEGIN
		cpp_dec_float_50 f_a = polynomial_eval(coef, start, N), f_b = polynomial_eval(coef, end, N);
	if (near_zero_2(f_a) || near_zero_2(f_b) || (f_a * f_b) < 0) {
		flag = true;
	}
	API_END

		return flag;
}

std::vector<cpp_dec_float_50> get_mon_interval_2(cpp_dec_float_50 coef[], cpp_dec_float_50 start, cpp_dec_float_50 end, int N) {
	std::vector<cpp_dec_float_50> intervals;
	intervals.push_back(start);

	// 递归将导数的零点放入
	cpp_dec_float_50 df[10];
	polynomial_derivative(df, coef, N);
	N = N - 1;
	std::vector<cpp_dec_float_50> df_roots = getroots_2(df, start, end, N);
	intervals.insert(intervals.end(), df_roots.begin(), df_roots.end());

	intervals.push_back(end);
	return intervals;
}


cpp_dec_float_50 getroot_2(cpp_dec_float_50 coef[], cpp_dec_float_50 start, cpp_dec_float_50 end, int N) {
	cpp_dec_float_50 root;

	//API_BEGIN
	if (start > end) {
		std::cout << "<<<<<<<<<<<<<<<<<<<<<<<<<<<<< Internal ERROR <<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl;
		//sys_error(API_FAILED);
		return root;
	}

	// 根在区间端点，直接返回
	if (near_zero_2(polynomial_eval(coef, start, N))) {
		return root = start;
	}
	if (near_zero_2(polynomial_eval(coef, end, N))) {
		return root = end;
	}
	//API_END

	// 根在区间内部
	API_BEGIN
		// 获取多项式的阶数

		int order = N;

	// 一阶多项式ax+b=0，根root = -b/a
	if (order == 1) {
		root = -coef[0] / coef[1];
	}

	// 二阶多项式ax^2+bx+c=0
	else if (order == 2) {
		//std::vector twice_roots = twice_order_roots(f);
		for (cpp_dec_float_50 x : twice_order_roots_2(coef, N)) {
			if (start <= x && x <= end) {
				root = x;
				break;
			}
		}
	}

	// 高阶多项式
	else {
		// 迭代初始点设置为中点
		cpp_dec_float_50 x = (start + end) / 2;
		cpp_dec_float_50 x_next;

		// 使用牛顿迭代法求根
		cpp_dec_float_50 df[10];
		cpp_dec_float_50 ddf[10];
		polynomial_derivative(df, coef, N);
		polynomial_derivative(ddf, df, N);
		for (int i = 0; i < 2000; i++) {
			cpp_dec_float_50 f_x = polynomial_eval(coef, x, N);
			cpp_dec_float_50 df_x = polynomial_eval(df, x, N);
			cpp_dec_float_50 ddf_x = polynomial_eval(ddf, x, N);

			// 函数值足够小视为零点
			/*if (near_zero_2(f_x)) {
				root = x;
				break;
			}*/
			//if (fabs(f_x) > 0.01) {
			//	x_next = newton_iter(f_x, df_x, x);
			//}
			//else {
			//	// 三阶改进牛顿法
			//	//x_next = three_order_newton(f_x, df_x, ddf_x, x, 15.0);

			//	//四阶牛顿法super halley
			//	x_next = four_order_newton_halley(f_x, df_x, ddf_x, ddf_x_2, x);

			//	//四阶牛顿法Jarratt
			//	//x_next = four_order_newton_Jarratt(f_x, df_x, df_y, x);
			//}
			x_next = newton_iter_2(f_x, df_x, x);
			//x_next = high_order_newton(f, x);
			//x_next = four_order_newton_halley(f_x, df_x, ddf_x, ddf_x_2, x);
			// 步长过小，迭代收敛
			if (near_zero_2(x_next - x)) {
				root = x_next;
				break;
			}
			//缩小区间
			cpp_dec_float_50 f_xmin = polynomial_eval(coef, start, N);
			cpp_dec_float_50 f_xmax = polynomial_eval(coef, end, N);
			cpp_dec_float_50 f_xnext = polynomial_eval(coef, x_next, N);
			if (f_xnext * f_xmin > 0) {
				start = x_next;
			}
			if (f_xnext * f_xmax > 0) {
				end = x_next;
			}
			// 超出求根区间，使用牛顿切线法拉回
			if (x_next < start || x_next > end) {
				cpp_dec_float_50 f_a = polynomial_eval(coef, start, N);
				cpp_dec_float_50 f_x = polynomial_eval(coef, x, N);
				cpp_dec_float_50 f_b = polynomial_eval(coef, end, N);
				// 根在迭代点左侧
				if ((f_a * f_x) < 0) {
					x_next = start - f_a * (x - start) / (f_x - f_a);
				}
				// 根在迭代点右侧
				else {
					x_next = end - f_b * (x - end) / (f_x - f_b);
				}
			}
			x = x_next;

		}
		root = x;		// 否则陷入循环的时候会返回未初始化的root值
	}
	API_END
		return root;
}


std::vector<cpp_dec_float_50> getroots_2(cpp_dec_float_50 coef[], cpp_dec_float_50 start, cpp_dec_float_50 end, int N) {
	std::vector<cpp_dec_float_50> roots;

	API_BEGIN
		// 获取多项式的阶数

		int order = N;
	//std::cout << "现在是" << N << "次方程" << std::endl;
	// 一阶多项式ax+b=0，根root = -b/a
	if (order == 1) {

		cpp_dec_float_50 root = -coef[0] / coef[1];
		if (start <= root && root <= end) {
			roots.push_back(root);
		}
	}

	// 二阶多项式ax^2+bx+c=0
	else if (order == 2) {
		for (cpp_dec_float_50 x : twice_order_roots_2(coef, N)) {
			if (start <= x && x <= end) {
				roots.push_back(x);
			}
		}
	}

	// 高阶多项式
	else {
		// 求单调区间
		std::vector<cpp_dec_float_50> intervals = get_mon_interval_2(coef, start, end, N);

		//std::cout << "---------------------------------已分区间---------------------------------" << std::endl;

		/*for (int i = 0; i < intervals.size() - 1; i++) {
			std::cout << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10) << "(" << intervals[i] << "," << intervals[i + 1] << ")" << std::endl;
		}*/



		// 对每个区间求根
		for (int i = 0; i < intervals.size() - 1; i++) {
			if (hasroot_2(coef, intervals[i], intervals[i + 1], N)) {
				roots.push_back(getroot_2(coef, intervals[i], intervals[i + 1], N));
			}
		}
	}

	API_END
		return roots;
}





// ************************************************************
outcome answer_api_nroots_of_law(law* input_law, double start, double end, int* size, double** answer) {
	API_BEGIN
		if (api_checking_on) {
			check_law(input_law);
		}
	// 转换为多项式
	//law_polynomial* f = input_law->polynomial(input_law);
	// 输出信息
	//f->print();
	//std::cout << "degree:\t" << f->degree() << std::endl;
	//std::cout << "coefficient:\t0:" << f->coefficent(0) << "\t1:" << f->coefficent(1) << "\t2:" << f->coefficent(2) << std::endl;
	// 表达式
	//std::cout << "in() = \t" << f->in()->string() << std::endl;
	//std::cout << "factor() = \t" << f->factor()->string() << std::endl;
	//std::cout << "expand() = \t" << f->expand()->string() << std::endl;

	// 输出表达式
	//std::cout << "f(X) = \t" << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10) << input_law->string() << std::endl;
	// 输出求导表达式
	//std::cout << "df(X) =\t" << input_law->derivative(0)->string() << std::endl;
	// 输出表达式X=1的值
	//std::cout << "f(0) = \t" << input_law->eval(0.0) << std::endl;

	//std::cout << "roots: " << std::endl;
	std::vector<double> roots = getroots(input_law, -100.0, 100.0);

	//for (double x : roots) {
		//std::cout << x << "\t";
		//double f_test = input_law->eval(x);
		//std::cout << "f(x)= " << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10) << f_test << std::endl;
	//}
	//std::cout << std::endl;

	//double f_real = input_law->eval(-1.992452178949137);
	//std::cout << "f_real= " << f_real << std::endl;
	// 除去重复值，返回结果
	std::sort(roots.begin(), roots.end());					// 排序
	auto last = std::unique(roots.begin(), roots.end());	// 将相邻且重复的元素踢到尾部，返回第一个重复元素的指针
	roots.erase(last, roots.end());							// 删去尾部的重复元素
	*size = roots.size();
	*answer = ACIS_NEW double[roots.size()];
	memcpy(*answer, &roots[0], roots.size() * sizeof(double));

	API_END
		return result;
}

outcome answer_api_nroots_of_coef(cpp_dec_float_50 coef[], cpp_dec_float_50 start, cpp_dec_float_50 end, int* size, cpp_dec_float_50** answer, int N) {
	API_BEGIN

		//std::cout << "roots: " << std::endl;
		std::vector<cpp_dec_float_50> roots = getroots_2(coef, -100.0, 100.0, N);

	//for (cpp_dec_float_50 x : roots) {
		//std::cout << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10) << x << "\t";
		//cpp_dec_float_50 f_test = polynomial_eval(coef, x, N);
		//std::cout << "f(x)= " << f_test << std::endl;
	//}
	//std::cout << std::endl;

	// 除去重复值，返回结果
	std::sort(roots.begin(), roots.end());					// 排序
	//auto last = std::unique(roots.begin(), roots.end());	// 将相邻且重复的元素踢到尾部，返回第一个重复元素的指针
	auto last = std::unique(roots.begin(), roots.end(), PrecisionComparator(1e-15));
	roots.erase(last, roots.end());							// 删去尾部的重复元素
	/*for (int i = 0; i < roots.size(); i++) {
		std::cout << std::setprecision(std::numeric_limits<cpp_dec_float_50>::digits10) << "剔除后：" << roots[i] << "\t";
	}*/
	*size = roots.size();
	*answer = ACIS_NEW cpp_dec_float_50[roots.size()];
	memcpy(*answer, &roots[0], roots.size() * sizeof(cpp_dec_float_50));

	API_END
		return result;
}
